% 用来将预测的结果可视化

%% 初始化环境
clc
clear
close all
landcolor = [0.98 0.97 0.97];
lakescolor = [0.71 0.84 0.92];
 
%% 从经纬度出发
[~,~,result]=xlsread('result.xlsx','Sheet1');
[m,n] = size(result);
 
%% 将含有经纬度的数据点画在地图上
mcolor = [249 222 239; 255 86 0; 236 165 201; 47 70 150];
mcolor = mcolor ./255;
 
%% 开始绘图
figure
 
%画出地图的经纬度范围
ax = worldmap('World');
ax = worldmap([40, 80], [220, 280]);
 
%画出美国大陆地图
land = shaperead('landareas', 'UseGeoCoords', true);
geoshow(ax, land, 'FaceColor', landcolor)
 
%画出美国湖泊地图
lakes = shaperead('worldlakes', 'UseGeoCoords', true);
geoshow(lakes, 'FaceColor', lakescolor)
 
%画出美国各个州，根据人口密度确定不同颜色
states = shaperead('usastatelo', 'UseGeoCoords', true);
scolor = [];
for i = 1:size(states, 1)
     if states(i).PopDens2000 >= 150    %获取人口密度
        scolor = [scolor; mcolor(4, :)];
     elseif states(i).PopDens2000 >= 50
        scolor = [scolor; mcolor(3, :)];
     else
        scolor = [scolor; mcolor(1, :)];
     end
end
index = find(scolor(:, 1) >= 1);
for i = 1:size(index, 1)
     scolor(index(i), :) = [1 0 0];
end
faceColors = makesymbolspec('Polygon',...
    {'INDEX', [1 numel(states)],  'FaceColor', ...
    scolor}); 
geoshow(ax, states, 'DisplayType', 'polygon', ...
 'SymbolSpec', faceColors)
 
%画出美国内陆河流
rivers = shaperead('worldrivers', 'UseGeoCoords', true);
geoshow(rivers, 'Color', 'blue')

%%
colorset=[240.0 255.0 11.0 ;...
233.38235294117646 247.73529411764707 13.5 ;...
226.76470588235293 240.47058823529412 16.0 ;...
220.14705882352942 233.2058823529412 18.5 ;...
213.52941176470588 225.94117647058823 21.0 ;...
206.91176470588235 218.6764705882353 23.5 ;...
200.29411764705884 211.41176470588235 26.0 ;...
193.6764705882353 204.14705882352942 28.5 ;...
187.05882352941177 196.88235294117646 31.0 ;...
180.44117647058823 189.61764705882354 33.5 ;...
173.8235294117647 182.3529411764706 36.0 ;...
167.2058823529412 175.08823529411765 38.5 ;...
160.58823529411765 167.8235294117647 41.0 ;...
153.97058823529412 160.55882352941177 43.5 ;...
147.3529411764706 153.29411764705884 46.0 ;...
140.73529411764707 146.02941176470588 48.5 ;...
134.11764705882354 138.76470588235293 51.0 ;...
127.5 131.5 53.5 ;...
120.88235294117648 124.23529411764707 56.0 ;...
114.26470588235296 116.97058823529412 58.5 ;...
107.64705882352942 109.70588235294119 61.0 ;...
101.02941176470588 102.44117647058823 63.5 ;...
94.41176470588238 95.1764705882353 66.0 ;...
87.79411764705884 87.91176470588235 68.5 ;...
81.1764705882353 80.64705882352942 71.0 ;...
74.55882352941177 73.38235294117646 73.5 ;...
67.94117647058823 66.11764705882354 76.0 ;...
61.323529411764724 58.85294117647058 78.5 ;...
54.70588235294119 51.58823529411765 81.0 ;...
48.08823529411765 44.323529411764696 83.5 ;...
41.470588235294144 37.05882352941177 86.0 ;...
34.85294117647061 29.794117647058812 88.5 ;...
28.235294117647072 22.529411764705884 91.0 ;...
21.617647058823536 15.264705882352956 93.5 ;...
15.0 8.0 96.0 ];

%%
for i = 2:36
    pointcolor=[colorset(i-1,1),colorset(i-1,2),colorset(i-1,3)];
    point(1).Geometry = 'Point';
    point(1).Lon = result{i,4};
    point(1).Lat = result{i,3};
    point(1).Name = 'primary';
    geoshow(point,'MarkerEdgeColor', pointcolor/255, 'Marker', '.', 'MarkerSize',10+i/3-2 )
end

%矫正图形
set(gcf, 'Position', [100, 100, 800, 400])